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ABSTRACT 


Matched-Field Processing (MFP) and Matched-Mode Processing (MMP) 
are two popular techniques for passively localizing an underwater acoustic 
emitter in range and depth. One major drawback of these techniques has been 
their sensitivity to uncertainty concerning the acoustic environment. Several 
methods of addressing this phenomenon have been proposed in the literature, 
with varying degrees of success. Achieving high-quality location estimates 
remains a problem except in simple range-independent experiments or 
numerical simulations. In this study, we demonstrate an approach for robust, 
accurate emitter localization in a highly range-dependent real environment 
using MMP. The main factors contributing to successful localization are: 1) use 
of the high-resolution Multiple Signal Classification (MUSIC) algorithm, which 
performs well even when only a few robust modes can be obtained by mode 
filtering; and 2) use of an acoustic propagation model incorporating mode 
coupling, which is able to generate accurate replica fields in a strongly range- 
dependent environment. A secondary objective of the study was to 
demonstrate the application of higher-order statistical estimation techniques 
to reduce noise effects. Our results indicate that these techniques show 
unacceptable sensitivity to noise- and model-induced estimation errors and 
require further refinement before they will be useful in the underwater acoustic 


localization problem. 
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I. INTRODUCTION 


Throughout the history of anti-submarine warfare, there has been 
great interest in the use of passive acoustic measurements to localize 
submarines. Numerous methods of estimating the Direction-Of-Arrival 
(DOA) of acoustic emissions from a target of interest have been developed 
over the years for sonar applications. However, these techniques are 
inherently incapable of directly determining the range and depth of a target 
(emitter) although there are indirect means of determining range by 
observing target DOA as a function of time. Because knowledge of range and 
depth is so vital in military applications, there has been considerable interest 
in developing techniques for direct determination of these parameters. One 
such technique, which has attracted considerable attention in recent years, is 
a generalization of DOA estimation known as Matched-Field Processing 
(MFP), along with a variation on MFP known as Matched-Mode Processing 
(MMP). 


Because of the similarity between DOA estimation and MFP/MMP, 
many of the techniques used in DOA estimation may be generalized for use in 
MFP/MMP. Two of the most popular DOA estimation methods—Bartlett and 
Minimum-Variance—have been studied extensively in the context of MFP 
(although the Minimum-Variance method has not been addressed in MMP) 
[Refs. 1, 2, 3, 4, 5, 6]. The MUSIC method has received extensive coverage in 
the DOA estimation literature, but relatively scant attention in the MFP 
literature [Refs. 7, 8, 9] and no attention in the MMP literature, possibly 
because of a perception that it would not perform well in realistic underwater 


acoustic environments. 


One of the most noteworthy results from MFP/MMP studies to date is 


the great sensitivity of the location estimates to uncertain knowledge of the 


acoustic environment. Numerous researchers [Refs. 10, 11, 12, 13, 14, 15] 


have studied this phenomenon and proposed various methods for addressing 


it, with varying degrees of success. The problem remains an open issue. 


Because of this sensitivity, and also because of the limited availability of real 


data sets, most MFP/MMP research to date has involved either simple range- 


independent experiments or numerical simulations. 


The objectives of this dissertation are to: 


Demonstrate the application of MFP/MMP to experimental data 
obtained in a strongly range-dependent environment during the 


1992 Barents Sea Polar Front Experiment [Refs. 16, 17, 18, 19]; 


Demonstrate that a coupled-mode propagation model is able to 
model the acoustic fields used in MFP/MMP with sufficient 


accuracy for high-quality localization estimates; 


Demonstrate that the high resolution of the MUSIC algorithm in 
combination with MMP is able to produce accurate location 


estimates in a realistic environment; and 


Demonstrate the application of higher-order statistics to the 
MFP/MMP problem. Although such methods have received some 
attention in the DOA estimation literature [Refs. 20, 21, 22], they 
have not yet been applied to MFP/MMP. 


Chapter II gives an overview of pertinent background material, 


including DOA estimation algorithms, modeling of underwater sound 
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propagation, and MFP/MMP theory. In that chapter, we: describe three of the 
most popular DOA estimation algorithms—the Bartlett, Minimum-Variance, 
and MUSIC methods; review the modeling of the acoustic field via 
decomposition into normal modes for both range-dependent and range- 
independent environments; discuss the application of higher-order statistics 
to DOA estimation; and derive the extension of DOA estimation techniques to 
MFP/MMP. Except for the portion regarding application of the MUSIC 
algorithm and higher-order statistics to MMP, this chapter contains no 
original material. Chapter III gives a brief overview of those features of the 
Barents Sea Polar Front Experiment which are relevant to this dissertation, 
including the physical characteristics of the channel and a description of the 
emitter and receiver. Chapter IV provides additional information concerning 
the data analysis procedure; it also presents and interprets the results of the 
analysis. Chapter V gives the conclusions reached from the research and 


proposes areas for further investigation. 





II. BACKGROUND 


This chapter provides the framework for our analytical approach. It 
includes overviews of the following topics: DOA estimation fundamentals, 
including an extensive description of the MUSIC algorithm; application of 
higher-order statistics to DOA estimation; normal mode modeling of the 


acoustic field; and theoretical foundations of MFP/MMP. 


A. DOA ESTIMATION 


One of the most fundamental parameters of interest in military 
applications is the Direction-of-Arrival (DOA) corresponding to a target of 
interest. This parameter is one of the primary outputs of virtually all military 
radar and sonar systems. DOA may be expressed in terms of azimuth 
(bearing) and/or elevation. Over the years, numerous technigues have been 
developed for DOA estimation (for a good overview, see [Ref. 23] and 
references therein); the most important of these techniques are discussed 
briefly in this chapter. As we will see later, these techniques may be 
generalized in a straightforward manner for use in Matched-Field Processing. 

1. Signal Model 

As is generally the case in signal processing algorithms, DOA 
estimation techniques make certain assumptions about the signals being 
processed. They assume, in particular, that the sound pressure field in the 
underwater acoustic channel may be expressed as a plane wave (i.e., that the 
surfaces of constant phase are planar). As we will see later, this assumption 
is not true in general, but is useful in many cases of practical interest. We 
will also assume that the signal is temporally narrowband with center 


frequency @; i.e., that amplitude and phase modulation do not introduce 


appreciable amplitude and phase changes over the physical extent (length) of 
the receiving array. If a signal does not satisfy this latter condition, the 
signal may be decomposed via Fourier techniques into narrowband 
components which do satisfy the condition. For the sake of generality, during 
most of the background discussions, we will allow the number of emitters to 
be arbitrary, even though the presence of a single emitter will be assumed 
during all actual data analysis. 

For mathematical convenience, we will conduct our analysis using the 
complex envelope representation of the signal rather than the real (physical) 
signal itself. Thus, if s¢) is the real signal, the corresponding analytic signal 


or pre-envelope [Ref. 24] is given by 
s(t) = s(t) + js(t), 
where s(t) is the Hilbert transform of s(t), 


S 
iz 


s(t) = n Was, 
The pre-envelope of a real bandpass signal at center frequency @ may be 
determined as follows [Ref. 25]: 
e Multiply the real signal s(t) by the complex carrier exp(j œt); 


e Pass the resulting signal through a low-pass filter to remove the 


component at twice the carrier frequency; 
e Multiply the resulting baseband signal by the complex carrier. 


All signals appearing in the sequel will be the pre-envelopes of real 


narrowband signals unless otherwise indicated. The narrowband assumption 


mentioned earlier is equivalent to requiring the complex pre-envelope of the 


received signal to be of the form 
s(t) 2 S, exp(Jax), 


where the (complex) amplitude S; is a slowly varying function of time, i.e., 


S,(t) = s[-5). 


U 


where / 1s the length of the array and v is the speed of sound. 

The concept of an array response vector is one that arises often in DOA 
estimation and MFP. To illustrate this concept, we consider an arbitrary 
receiving array of N elements. We may represent the signal as an N- 
dimensional, time-varying, complex vector whose components are the signals 
at the individual array elements. For convenience, we will assume that only 
the azimuth @ of the target is of interest, although extension to include 
dependence on elevation is straightforward. In DOA estimation, the array 
response vector a(0) is defined to be the unit-normalized (i.e., length of 
a(@)=1), noise-free, pressure field vector expected (based on the modeling 
assumptions) at the receive array given that an emitter 1s at the angle 0. The 
set of a vectors for all possible values of 0 is known as the array manifold. 
More generally, a could be a function of parameters other than DOA as well. 
For the simple case involving a single emitter, a receiving array of N 
identical, omnidirectional elements in an arbitrary geometry and a 
narrowband, plane-wave signal with center frequency @, a(6) may be 


expressed as 


where 7, (a function of 0) represents the time delay seen by the incoming 
plane wave between sensor z and sensor 0 and superscript 7' denotes (non- 
conjugate) matrix transpose. In the general case where d emitters are 


present, the signal model used in DOA estimation is 
p=As+n, (1) 


where p=[p,(t) p(t) -~ p. (t)| and n=[n,(t) n,(t) -- ny (t)]- are 
vectors whose elements are the received pressure and noise, respectively, at 
each element of the array; s=|s,(¢) s(t) = ss] is the vector of signals 
produced by the d emitters (s,(t)2S;,exp(jot), with S; a complex amplitude, 
because of the narrowband assumption); A= la(6,) a(0,) -- a(0, )| is a 
matrix whose columns are the array response vectors corresponding to the 
DOAs of the d emitters; and ¢ is time. Assuming that the signal and noise are 
uncorrelated, the signal-plus-noise covariance matrix R, is then given by 
R, - F[pp" | 
- ARA" « R,, (2) 


where R, = E|ss" | and R,, = E|nn" | are the signal and noise covariance 
matrices, respectively, superscript H denotes Hermitian (conjugate) 
transpose and E denotes statistical expectation. An estimate of R, derived 
from the measured data is the fundamental quantity used in virtually all 
DOA estimation algorithms studied to date. 

2. Algorithms 

Numerous signal processing algorithms have been studied in the 
context of DOA estimation. Many are fairly obvious generalizations of 


techniques used for estimating the spectra of time series. 


a. Weighted-sum beamforming 

Weighted-sum beamforming (see, for example, [Ref. 26]) is the 
most commonly used technique (in practice) for DOA estimation, and is used 
in virtually all modern military radar and sonar systems. This technique is 
analogous to the Finite Impulse Response (FIR) filters used in time-series 
analysis. The output b(t) of the beamformer is simply a linear combination of 


the signals received by the elements of the array, 
b(t) = w” p(t), 


where w is the vector of weights. The weight vector w is chosen to satisfy a 
statistical constraint which is appropriate for the given situation. The 


expected value B of the output power of the beamformer is given by 


B= Elle] -w'R,w. (3) 


A particular value for w generally produces high gain only in a single look 
direction, i.e., for a single value of 0. In practice, multiple look directions are 
of interest, so that multiple w vectors are required (this approach is 
analogous to the use of a "bank" of matched filters in, for example, active 
sonar). Thus, in general, both w and B are functions of 0. The value of 0 for 
which B is maximized is taken as the estimated DOA of the emitter; this 
value may be obtained by conventional one-dimensional search techniques. 
The Bartlett beamformer is a special case of weighted-sum 
beamforming in which the weight vectors w are simply the array response 
vectors a for all look directions of interest. When the noise is spatially 
homogeneous, it can be easily shown that the output of this beamformer has 
the highest possible Signal-to-Noise Ratio (SNR) of any weighted-sum 


beamformer. This beamformer is analogous to the periodogram [Ref. 27] used 
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in time-series analysis. Since for the Bartlett beamformer w=a(@), the 


expected value of the output power of the beamformer is 
B(@) = a"(@)R,a(@) (4) 


and the DOA estimate is given by 


^ 


0 — arg max a" (0)R,a(0). 

The Minimum-Variance method (MVM), also (somewhat 
misleadingly) known as the “Maximum Likelihood” (ML) beamformer [Refs. 
28, 29] selects the beamformer weights to minimize the beamformer output 
power, subject to the constraint of unity gain (r.e., zero distortion) in the 
desired look direction 0. Formally, w 1s chosen to 

minimize w^ Row 


subject to w^a(0) — 1, 


The effect of this minimization is to produce the lowest array response at 
directions that have the strongest signals. Thus, this method is useful in 
reducing the effects of directional noise on beamformer performance. The 
required weights may be easily shown (e.g., [Ref. 23]) (using the method of 
Lagrange multipliers) to be 

R,a(0) 
a”(O)R,a(6) 


w(6)- 
Substituting into (3) and simplifying gives a beamformer output power of 

B(0) - EE (5) 
For a comparison of the performance of the Bartlett and MVM techniques, 


see Lacoss [Ref. 29]. 
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b. MUSIC 

The MUSIC (for MUltiple SIgnal Classification) method [Refs. 
30, 31, 32] 1s one of the first and probably the most widely studied of a class 
of techniques which address the DOA estimation problem from a geometric 
perspective; such methods offer potentially large improvements in resolution 
with respect to beamforming [Ref. 23]. Because MUSIC is central to the 
investigations of this dissertation, a full discussion of it is provided in the 
sequel. This discussion will attempt to develop an intuitive understanding of 
the algorithm rather than to provide a rigorous proof of its methods. 

Signal Subspace. Figure 1 is a geometric illustration of the 
behavior of the received signal vector p(t) due to d emitters at various DOAs. 
The coordinate axes represent the responses of the N sensors. The 
components of the vector p are the outputs of the individual sensors and are 
functions of time. As time progresses, the tip of p then sweeps out a curve as 
shown in the figure. Due to obvious limitations, our illustrations can show at 
most N=3 and will show sensor outputs as real. These limitations will not 


adversely affect the illustration of key concepts. 


Sensor 3 


Sensor 2 





Sensor 1 
Figure 1: Behavior of signal vector 


ET 


Now consider an example involving a single emitter (d=1; Figure 
2). For this case, as time progresses, each component of p is multiplied by the 
same time-dependent phase factor. The (complex) magnitude of p changes, 
but its direction does not. Therefore, the “curve” swept out by p in this case is 


a straight line passing through the origin. 


Sensor 3 


Sensor 2 





Sensor 1 
Figure 2: Signal Subspace (d=1) 
Next consider an example involving two emitters (d=2) at angles 
0, and 6, (Figure 3). For this case, p is the vector sum of contributions from 


the individual signals; specifically, 


p = s,(t)a(0,) + s,(¢)a(@,). (6) 

These contributions s,(/)a(0; are vectors whose magnitude 

varies with time, but whose direction is fixed by a(@,). Since the two vectors 
are multiplied by different time-varying phase factors s,(f)-S,exp(J%té) 
(provided that the slowly varying complex amplitudes S, are not perfectly 
correlated), their sum p sweeps out a curve which is confined to a plane 


passing through the origin. 
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In general, for d emitters, the tip of the received signal vector p 
sweeps out a curve which is confined to a d-dimensional subspace of (”. This 
subspace is known as the signal subspace. Intuitively, then, the number of 
emitters could be determined by measuring the dimension of the signal 


subspace. We will show later how this can be done. 


Sensor 2 





Sensor 1 
Figure 3: Signal Subspace (d-2) 

Array Manifold. Figure 4 illustrates an array manifold, 1.e., 
the locus of the array response vectors a(0) for all possible values of 0. Note 
that, by definition, a is independent of time, so that time does not participate 
in this illustration. In practice, it is not necessary to have an analytical 
expression for a(0); we can instead determine it experimentally at a finite 
number of points, store the results, and recall them when desired 


(interpolating when necessary). 
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It may happen that the array manifold “runs over itself” (i.e., the 
mapping a(@) from the interval 0 to 2z into €" is not one-to-one). In such a 
case, the array is said to have an ambiguity. Such an ambiguity is not the 
only kind possible. For example, when the a vectors corresponding to 3 
different DOAs lie in a single plane, the array also has an ambiguity (for 
reasons which may not be obvious at this point). In general, when the a 
vectors corresponding to n+1 different DOAs lie in a subspace of dimension n 


or less, the array is said to have a rank n ambiguity. 


Sensor 3 


V 


Sensor 2 


Sensor 1 
Figure 4: Array manifold 

DOA Determination. For the case of one emitter, it is clear 
that the array response a corresponding to the emitter DOA lies in the (1- 
dimensional) signal subspace illustrated in Figure 2. In order to determine 
the DOA, we would observe the behavior of the signal vector p (if no noise 
were present) to determine the signal subspace and find the single unit vector 
that spans the subspace. This unit vector represents a point on the array 


manifold corresponding to the angle of the emitter. We would then invert the 
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mapping a(0) (which is one-to-one unless the array has an ambiguity) to 
determine the DOA 8. 

Now consider the case of two emitters and recall that the signal 
vector p is the sum of (vector) contributions from the individual emitters (6). 
At any time £, these contributions are scalar multiples of the a vectors 
corresponding to the two emitter DOAs. Thus, the signal subspace is spanned 
by two vectors a(@,) and a(@,) which correspond to the points where the array 
manifold intersects the signal subspace (Figure 5). If the array has no 
ambiguities, there is no third a vector that lies in the subspace. Thus, just as 
in the case of one emitter, we can invert the a(@) mapping provided by the 


array manifold to determine the DOAs (Figure 5). 


Sensor 3 





Sensor 2 


Figure 5: DOA Determination for 2 emitters 


In general, we observe the data vector p, determine the signal 
subspace, find its intersection (d different a vectors) with the array manifold, 
and invert the mapping a(@) to estimate the DOAs of the emitters. It is 


apparent that, in general, this method requires the number of emitters d to 
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be fewer than the number of array elements N (if d2N, the entire array 
manifold may lie within the signal subspace). Clearly, if the array manifold is 
such that the a vectors corresponding to the emitter DOAs are not linearly 
independent (array ambiguity), the signal subspace dimension 1s less than 
the number of emitters. In such a case, the method described above fails to 
give the correct number of emitters and fails to identify one or more of the 
emitter DOAs. Such ambiguities can often be avoided by proper array design 
(although sometimes physical constraints, such as in line arrays, preclude 
such design). 

Multipath. It is possible for the emissions from a single emitter 
to arrive by two or more different paths (e.g., as a result of reflection from the 
water surface, in the sonar case). To illustrate the effect of such a situation, 
we consider the case of an array receiving signals from one emitter via two 
different DOAs. In such a case, the two signals will exhibit perfect temporal 
correlation (at least in theory). Recall that it is the independent variation 
with time of the array output due to individual signals from different DOAs 
that causes the signal vector to sweep out a two-dimensional subspace 
(Figure 3). In the multipath case, however, the array no longer responds 
independently to signals received from the two DOAs: the signal subspace is 
one-dimensional. In addition, note that this signal subspace is not spanned by 
any combination of vectors in the array manifold. Thus, the geometric 
approach is incapable of dealing with perfectly coherent multipath. 
Fortunately, such multipath is rarely, if ever, observed in real life. However, 
the performance of most geometrically-based algorithms degrades rapidly as 
the correlation between emitters (or multipath arrivals) approaches 1. 

Noise. The above discussion assumes that no noise is present. 


With noise present, the signal-plus-noise vector p sweeps out a curve that no 
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longer lies in the signal subspace. Thus, the simple technigues described are 
not valid if noise 1s present. 

Essentially, the technique described above for noise-free 
conditions must be modified as follows. For simplicity, we will assume that 


the noise 1s spatially 1sotropic and uncorrelated, so that (2) becomes 


H 
R, - AR,A" *c'I,, 
where I, is the NXN identity matrix and the covariance matrices R, and R, 


represent theoretical, not estimated statistics. We now consider the following 


eigenvalue problem: 


Since R, is Hermitian, the eigenvalues are real and the associated 


eigenvectors may be selected to be orthonormal [Ref. 33]. Then, 
0 = det[R,, - Al, | 
= det[AR,A” - (4-07) | (7) 


Per the definition of the Nxd matrix A, the columns of A are the a vectors 
corresponding to the d emitters. As mentioned above, if the array has no 
ambiguities, these a vectors are linearly independent and therefore A has full 
(column) rank d. The elements F.J) of R, have the form r,S,S*, where r, is a 
correlation coefficient between the ith and jth signals and S, represents the 
(complex) signal amplitude (not a function of time in this case) of the ith 
signal. Thus, the dxd matrix R, has full rank d unless one or more pairs of 
emitters are perfectly correlated (in which case r;-1 for some i, j and 


therefore the ith and jth columns are linearly dependent). Thus, assuming no 


array ambiguities and less than perfectly correlated emitters, the first term 
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in equation (7) is an NXN matrix of rank d. That term therefore has d non- 
zero eigenvalues v, and N-d zero eigenvalues. Since the eigenvectors of this 
term are also eigenvectors of the second term (identity matrix), the 
eigenvalues of the sum of the two terms in equation (7) are the sums of the 


respective eigenvalues; t.e., the eigenvalues of R, are 


Uile lvi +0°,..., Vy +07,07,… 07). 
The eigendecomposition of R, may thus be expressed as 
R, - EAE" 
- E;N,E$ Fol, 
where E is the matrix of eigenvectors, A=diag(A,), E, contains the columns of 


E corresponding to non-zero eigenvalues of AR,A", and Ng-diag(v). Now 


consider the equality 


ARA" — ENE. (8) 
Obviously, both the dxd diagonal matrix Ng and the Nxd matrix Eg have rank 
d. Thus we see that the right hand side of (8) consists of a NXN matrix with 
rank d, each of whose columns is a linear combination of the d columns of Eg. 
Therefore, the N columns of the right hand side must span the same subspace 
as the d columns of Eg. Similarly, the columns of the left hand side of (8) 
must span the same subspace as the d columns of A. Consequently, the 
columns of E, must span the signal subspace (the subspace spanned by the 
columns of A). In the presence of noise, then, the signal subspace may be 
determined by performing an eigendecomposition. The N-d eigenvectors 
corresponding to the zero eigenvalues of AR,A” span what is referred to as 
the noise subspace (the orthogonal complement to the signal subspace). In the 


case where the noise is not isotropic and spatially uncorrelated, a generalized 
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eigendecomposition | [Ref. 23] must be used; such a decomposition requires 
that an estimate of the noise covariance be available. The existence of an 
estimate of the noise covariance will not be assumed in the data analysis to 
follow, so we will not address this generalized eigendecomposition further. 
Obviously, there is some price to be paid for not incorporating the structure of 
the noise covariance in our method; we expect this price to increase as signal- 
to-noise ratio decreases. 

In the real world, the covariance matrices are never known 
exactly, but must be estimated from observations. In the data analysis to be 


presented later, the following estimate of the covariance matrix will be used: 
> _ H 
R, = PP TL. 


where L is the number of observations and P is an NxL matrix whose 


columns are observations of the received signal p at successive times; t.e., 
P - [p(&) | p(5) 7: p(t)]. 


Because of estimation errors, the N-d noise eigenvalues of this estimated 
covariance matrix will not have exactly the same value o^, so that even the 
number of emitters d cannot be estimated with certainty. For the purpose of 
the present discussion, however, we assume that d 1s known. Estimates of the 
signal subspace may then be obtained via the familiar Linear Least Squares 
and Maximum Likelihood techniques [Ref. 23]. For the simple case of 
Gaussian noise, these techniques give the same result, but one which 1s not 
computationally feasible in most practical situations (since both involve a 
global minimization over a space with dimension equal to the number of 
emitters). The MUSIC algorithm discussed in the sequel arose out of the need 


to reduce computational complexity. 
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Estimating the Subspace Dimension. In an optimal 
estimator, the problem of estimating the subspace dimension cannot be 
decoupled from that of estimating the subspace itself. However, in the 
interest of reducing computational load, we can estimate the dimension 
separately. This feature of the MUSIC algorithm therefore causes it to be 
suboptimal. 

Estimating the Signal Parameters. As mentioned earlier, in 
the presence of noise, we can no longer depend on precise intersections 
between the signal subspace determined from the estimated covariance 
matrix R, and the array response vectors corresponding to the signal DOAs. 
To estimate these DOAs, we must therefore determine those a vectors which 
are "closest" (1n some sense) to the signal subspace. There are several 
possible methods for doing this; we consider only the simplest method 
(Conventional MUSIC) here. 

Recall from the earlier discussion that when the theoretical (i.e., 
not estimated) covariance matrix R, 1s used to determine the signal and noise 
subspaces, the array response vectors corresponding to the signal DOAs span 
the signal subspace and are orthogonal to the noise subspace. Thus, the 


squared length of the projection of an a vector onto the noise subspace, given 


by 


Length squared - a" (0)E, E& a(0), 
will be zero when @ is one of the emitter DOAs. However, because of 
estimation errors (as well as because of inaccuracies in the signal model used 
to determine a(0), when the estimated covariance matrix R, is used to 
determine the signal and noise subspaces (i.e., to determine E,), this quantity 


will generally not be zero for any a. We must therefore search the array 
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manifold for the set of d array response vectors which result in the lowest 
value for the guantity. An alternative method for accomplishing this is to 


define the function 


ue 
a" (Q)E,Exa(0)- 


P, (0) 5 (9) 


The estimates of the signal DOAs then correspond to the peaks of this 
function. Although this function is similar in some respects to beamformer 
output power functions such as that in (5), the peak heights of this function 
do not necessarily provide any information about the power in their 


respective components. 


B. HIGHER-ORDER STATISTICS 


In recent years, there has been increasing interest in the use of higher- 
order (i.e., order greater than 2) statistics in signal processing. One reason for 
this interest is that the statistics known as cumulants are identically zero for 
Gaussian random variables, provided that the order of the cumulant is 
greater than 2. Intuitively, we would expect that, in situations where the 
noise is Gaussian but the signal is not, cumulant-based methods offer 
potentially large performance improvements over conventional methods 
based on second-order statistics, by removing the noise without affecting the 
signal. A detailed treatment of higher-order statistics is beyond the scope of 
this dissertation; for the purposes of this discussion, a brief consideration of 
the 4th-order cumulant will suffice. Broader treatments of the topic 
(including the material in the sequel) may be found in articles by Shiryaev 
[Ref. 34], Brillinger [Ref. 35], and Brillinger and Rosenblatt [Ref. 36]; tutorial 
articles by Nikias and Raghuveer [Ref. 37] and Mendel [Ref. 38]; and the 
recent textbook by Nikias and Petropulu [Ref. 39]. 
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The Joint characteristic function of a set lx,, x;,..., x,] of n real 


random variables 1s defined [Ref. 25] by 


®(@,, %,...,o,) - Eُ[exp[j(öıxı + @z, +:::+ o,z,)]), 


where E denotes statistical expectation, as before. The form of this function 
obtained by making the substitution s,=j@, is known as the moment- 
generating function; the moments of the x, can be obtained from the 
coefficients of the Taylor expansion of the moment-generating function about 
s =0. The second characteristic function WV of these same random variables is 
defined as W-lIn ®. The joint cumulants of order r=k,t+k,+~+k, of these 
random variables are defined [Ref. 25] as 


dT ER 


kı kz eee kn 
Ja, JO; Jo, PA -0„-. : OO, -0 


Crom IE (-J) (10) 
i.e., the coefficients of the Taylor expansion of V about 0,20. For the 4th-order 
cumulant of zero-mean real random variables x;, it can be shown [Refs. 34, 


39] that (10) reduces to 


Cum(x,, Xo, %3,%,) 2 Elx,x,x,x,] —- Ebxx,l -Elxx,l 


-Elx,x,] - Elxsx,] - Elxyx,] - Elxsxl. (11) 


It may further be shown [Ref. 25] that, if the x, are jointly Gaussian, the 4th 


order moment is given by 
EM SE |x. | E|x,x, | + E|xyx; | E|x,x, | - E[xyx, |- EXE (12) 


therefore the cumulant vanishes as claimed. For complex random variables, 


(11) takes the form [Ref. 40] 
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go xxx Elx x]: Elx] 


AE (18) 


where the third term is generally assumed to be zero due to the symmetry 
property between the real and imaginary parts of a stationary, bandpass, 


complex process [Ref. 41]. We will retain this term for generality. 


A spatial 4th-order cumulant matrix C, may be defined as follows: 


[P.(t)p;(t), p,(t)p;(t),---, pu(t)ps(t)]!, (14) 


e O 


where “*” denotes complex conjugate and the p,(f#) are sensor signals (complex 
pre-envelopes) at each of N sensors in an array. This definition differs from 
the one used by Nikias and Petropolu (Ref. 39], but was selected so that C, is 
Hermitian. By substituting (11) into (14), and using vector notation, we 


obtain 


C,- El(pep )P o p) }- Elpep JEl(p p') | 


“Elpp” Elp'p']-Elpp' EP PT], (15) 


where o denotes the element-by-element product of the vectors and 

superscript T indicates (non-conjugate) matrix transpose. The fourth term of 

(15) 1s generally assumed to be zero due to the symmetry property discussed 

in conjunction with (13) but will be retained for generality. Substituting the 

signal model of equation (1) (for the case of a single emitter) into (15) gives 
C- Ey (as » a's' (asc as) — Elas o a's' E (as oa s "| 


-Elass'a"| o Ela's'sa”) = Elassa”) o Ela'ss'a”) 
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provided that the signal and noise are independent (the single-emitter 
assumption is made here for notational convenience and will also be made 
during the actual data analysis in the seguel). Since a is deterministic, it may 
be factored out of the expectation operators, and we obtain, after some 


algebra 
C, - y(a»a' (asa) . (16) 
where 
y » Cum[s(t), s' (t), s(£), s' (t) 


is the kurtosis measure of the (single) emitter signal s(£) and a 1s the array 


response vector corresponding to the location of the emitter. 


The structure of C, 1s therefore identical to that of the covariance 
matrix defined in (2), except that: 1) the noise covariance vanishes due to the 
properties of the 4th-order cumulant; 2) the array response vector a 1s 
replaced by aea', which is real; and 3) C, has rank 1 (see (16)) due to the 
assumption of a single emitter and is real. This structure allows the use of a 


modified version of the previously discussed MUSIC algorithm as follows: 


e Form an estimate C, of the C, matrix from the measured data 


using (15), with expectation operators replaced by time averages; 


« Perform an eigendecomposition of C 4- Because it has rank 1, there 


will be (theoretically) a single non-zero eigenvalue. 


e Select the eigenvector corresponding to the largest eigenvalue; as 
shown earlier (8), this eigenvector must span the same subspace as 


the a vector corresponding to the actual emitter location. The 


24 


remaining eigenvectors span the noise subspace and form the 


columns of Ej. 


e Form the function 


1 


log = H c 
la(6)oa*(6)] E,EX[a(@)oa *(6)| 


(17) 
the estimate for the emitter location corresponds to the peak of this function. 
This method can be easily generalized to the multiple-emitter scenario, 
provided that the signal subspace dimension is selected to correspond to the 
number of emitters. However, our experimental work with real data does not 


require this generalization. 


Cumulant-based versions of the Bartlett and MVM methods also exist 
[Ref. 39], but these will not be used in the data analysis. It should also be 
noted that the form of the 4th-order cumulant appearing in (15) is a reduced 
form of that used by Porat and Friedlander [Ref. 20]. 


C. ACOUSTICS AND MODELING 

1. Helmholtz Equation 

Whereas in DOA estimation the received signals are assumed to be 
plane waves, in MFP the pressure field amplitude p(7,z) in an underwater 
channel (see Figure 6) due to a narrowband emitter with center frequency w 


is assumed to satisfy the Helmholtz equation [Ref. 42] 


1 3 | op(r,z)| o*p(r,z) p 
M emcee) |g ENTS) = 18 
E r 5 j EE t k(r,z) p(r,z) - 0, (18) 


where r is the range from the emitter to the observation point, z is the depth, 


and kzo/c is the wavenumber. Cylindrical coordinates are used in (18) 
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because of the symmetry which exists when changes in sound speed in the 
azimuthal direction are negligible. The seguel will present an overview of 
methods for solving (18) numerically which are pertinent to this dissertation. 


In this section, it will be assumed that no observation noise is present. 


ro 
r0 EE 
z-0 i. 
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Figure 6: Generic Underwater Channel 


2, Normal Mode Solution 
In a channel where all properties (sound speed, water depth, bottom 
type, etc.) are independent of horizontal range, it is well known [Ref. 42] that 
the pressure field p at range r (relative to an arbitrary origin; see Figure 6) 
sufficiently far away from the emitter and at depth z may be expanded in 
terms of normal modes as follows: 
Diner S za AE BE "7l (19) 
BART TR) 
where: 2, and r, are the emitter depth and range, respectively; r is taken to 
be greater than r,; A is a constant which depends on the emitter power; and 


the Z, and k, are the eigenfunctions and eigenvalues, respectively, of the 


ordinary differential equation 
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dZ w 
el SR. re (20) 


The boundary conditions for (20) depend on the acoustic properties of the 
surface (which is always assumed to be pressure-release; i.e., {z=0)=0) and 
the bottom. Attenuation due to sediment and water column is incorporated, 
as is customary, via a small imaginary part in the &,,. Although the sum in 
(19) is over values of m from 1 to infinity, all modes for which m is greater 
than some integer M are attenuated enough to be ignored at the ranges of 
interest; such modes include the strongly bottom-interacting and evanescent 
modes. In the discussions to follow, we will ignore these modes and 
incorporate only the lowest M modes in our normal mode expressions. It 
should be noted that, in general (i.e., for range near zero), the expansion (19) 
is only approximate, since it only accounts for the discrete spectrum of the 
modal solution to the Helmholtz equation (18). This fact does not present a 
problem, however, because the contribution to the pressure field from the 
continuous spectrum is negligible at the ranges of interest in MFP/MMP. 
Numerous methods exist for implementing (19) on a computer (see, for 
example, [Kef. 43]). As is characteristic of solutions of boundary value 
problems such as (20), the Z,, are orthogonal with respect to the density 
function p(z), l.e., 


a Û Z, (z)Z, (z)dz 7 v,ó(m - n), (21) 


where 


= x eae 
val, (2) )d 


This orthogonality will prove useful later in the discussion of Matched-Mode 
Processing. 

3. Adiabatic Approximation 

As mentioned above, the normal mode expression (19) is strictly valid 
only in range-independent environments; nevertheless, it can be modified 
slightly to apply to a limited class of range-dependent environments. 
Whenever range dependence exists (due to a sloping bottom or change in 
sound speed profile, for example), it is clear that the Z,, and &,, must be 
functions of range. If the range dependence is relatively weak, it can be 
assumed that a mode does not exchange energy with other modes, but merely 
adapts itself to local environmental conditions. This assumption is known as 
the “adiabatic” approximation. In such a case, (20) is solved at each of the 
ranges of interest. The resulting Z, and &, are then functions of range, so 


that the pressure field can be expressed as 


oo stra y-y Fîn EER. (zem joe jfk rar! (22) 
m=14/ k OE To) To 
4. Coupled Mode Model 
Due to the energy exchange between modes in a strongly range- 
dependent environment, the use of a normal-mode model in such an 
environment becomes more complex. For simplicity, we assume that 
horizontal refraction and azimuthal scattering are negligible (see [Ref. 44] for 
a fully three-dimensional treatment). As will be discussed in more detail 
later, this approximation is satisfactory for the acoustic environment 
addressed in the data analysis to follow [Ref. 19]. 


Mode coupling may be accounted for in a manner consistent with the 


notation used above by defining a range-dependent mode amplitude A,,(r) 
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(1.e., the complex constant A in (22) becomes a function of both range and 


mode number 7) [Ref. 45]. The pressure field is then given by 


u PO, TED, 
EE ES) E os OE) Zan eZ: forse jor- jÍ k,( oar (23) 


To 


In this expression, all coupling between modes is accounted for by the A, (i.e., 
each A, depends on the mode amplitudes and phases of all other modes along 
the propagation path). The Broadband Coupled Mode model developed by 
Chiu et al. [Ref. 19] has demonstrated high accuracy in predicting the modal 
structure observed in the Barents Sea Polar Front Experiment and will be 


used for all data analysis in the sequel. 


D. MATCHED-FIELD AND MATCHED-MODE PROCESSING 


A good overview of this topic may be found in the textbook by Tolstoy 
[Ref. 46]. Only background material pertinent to the later data analysis 
(along with some preliminary material) will be presented here. 
1. Motivation 
Several limitations in applying DOA estimation techniques to the 
underwater acoustic localization problem arise from the inability of the 
simple plane-wave signal model to describe the acoustic field adequately. In 
this section, we describe these limitations, which provide the motivation for 
study of MFP/MMP. 
a. Ability to determine target parameters 
Because of the assumption that the received signal is a plane 
wave, the only degree of freedom in the array response vector is DOA. 
Therefore, DOA estimation is inherently incapable of providing any other 


information. In particular, it gives no range information, which in military 
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applications is of vital interest (although it should be pointed out that it is 
possible to determine range by observing DOA information over time, 
provided that certain assumptions concerning target motion are valid; see 
[Ref. 47]). 

b. Estimation error 

Since the speed of sound in seawater is a function of position, 
the sound "rays" (1.e., paths normal to the surfaces of constant phase) are 
curved. Consequently, the DOA estimates will, in general, be different from 
the actual directions of the emitters. In many practical situations, the 
receiving array has no depth extent (e.g., a horizontal line array) and the 
ocean may be considered to be horizontally stratified (i.e., speed of 
propagation is a function primarily of depth). In such situations, the plane 
wave signal model is relatively accurate. Even in such cases, the DOA 
estimates will often be significantly in error due to reflection of sound energy 
from boundaries (surface and bottom). 

c. Loss of gain 

We have noted earlier that the Bartlett beamformer results in 
the highest output SNR, provided that the noise is spatially homogeneous. 
Essentially, this feature is due to the fact that the signals from different 
receive array elements add constructively (because the beamformer 
introduces phase shifts to account for the shape of the wavefront) while the 
noise components do not. In an underwater acoustic channel, however, the 
wavefronts may not be planar (particularly in the vertical direction), so that 
the signals from the array elements will no longer be added perfectly 
constructively (since the beamformer assumes the wavefront is planar when 
it actually is not). Consequently, the gain of the beamformer will be lower 


than when the incoming signal is a plane wave. 
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d. Inability to resolve multipath arrivals 

As mentioned earlier, many eigenstructure-based methods 
(including MUSIC) are unable to handle signals from multiple highly 
correlated emitters. Reflections from boundaries in the underwater acoustic 
problem can often result in the same signal arriving at the receiver from 
multiple directions. This situation is equivalent to the presence of multiple, 
highly correlated emitters; thus, MUSIC and similar methods break down in 
such a situation. 

2. Generalization of array processing algorithms 
a. Matched-Field Processing (MFP) 


For the case of a single signal, (1) becomes 
p=as(t)+n. (24) 


We can use (23) to express the pressure at a receiving array consisting of a 


set of N vertically-aligned hydrophones at depths {z,, Z,,...,Z,} 1n the form (24) 


if we identify 
PET) 
pel ze) (25) 
DOE) 
and 
aee EE. MASSE i ihrer] (26) 
e Nm 5) A 


where exp(jaX)-s(t) and Z, (r)-IZ, (rz), Z,(rzj,..., Z, (rz, . The array 
response vector a is now a function of emitter range r, and depth z, rather 


than azimuth, as in DOA estimation. The a(ry2,) are obtained by calculating 
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the pressure field seen at the receive array using an appropriate propagation 
model, for every possible (r,z,) combination of interest. In the previous 
discussions on DOA estimation, no assumptions were made concerning the 
form of a. Furthermore, the dependence of a on azimuth alone was for the 
purposes of illustration only, and is not required for the DOA estimation 
techniques to be valid. Therefore, the algorithms discussed above in the 
context of DOA estimation may be used in MFP as well [Refs. 3, 7]. It should 
be noted that MFP does not require that the acoustic field be expressed in 
terms of normal modes; the above analysis is valid for other types of 
propagation models as well. 

The well-known principle of acoustic reciprocity (see, e.g., [Ref. 
48]) is very useful in minimizing the computations required to generate the 
array manifold. This principle states that, under a set of reasonable 
assumptions, the acoustic pressure at a location (r,z) generated by a simple 
source at location (ry2,) is the same as the acoustic pressure at location (rz) 
generated by that same source at location (r,z). In MFP, the construction of 
the array manifold requires that the field at a known receiver location (r,z) be 
computed for every possible emitter location (7,,z,), whereas one run of a 
propagation model will generally produce the acoustic field at every possible 
receiver location due to an emitter at a known, fixed location. Thus, it may be 
seen that this principle of acoustic reciprocity allows construction of one 
component (corresponding to one element of the receive array) of the array 
manifold with a single run of a propagation model. The total number of runs 
required will be the same as the number of elements in the receive array. 
Clearly, this approach allows huge savings in computation compared to a 
“brute force” approach which performs one run of a propagation model for 


each possible source location. 
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The procedure for MFP may be outlined as follows: 


e Determine the pre-envelopes of the received (real) signals at each 
receive hydrophone. Use these pre-envelopes to generate an 
estimate (using time averages instead of statistical expectations) of 


the covariance matrix R, (or the 4th order cumulant matrix C,); 


e Using a suitable propagation model and invoking the principle of 
acoustic reciprocity, generate the array manifold vectors a(r,,z,) 
(which in the context of MFP/MMP are generally called replica 


fields) for values of (7,,z,) on a suitable grid; 


e Generate the functions in equations (4) (Bartlett), (5) (MVM), (9) 
(MUSIC), or (17) (cumulant MUSIC), as desired. These will now be 


functions of two variables (7,,z,) rather than one (9), i.e., 


B(r,, z,) » a” (n, z R a(n, zo) (Bartlett), (27) 


1 


e 
Uo» zo) a” (nZ) R aln, 2o) 


(MVM), (28) 


1 
D23 ES H 


————— m (MUSIC), 29 
a (nz E EZa(n 2) | ) e 


1 
a(n, z,)» a (ry, zy JE E,E5 |a(7, ze ne, )| 


The estimated emitter locations are the (ry2,) combinations for 


DM (Cumulant). (30) 


which these functions are maximized. In MFP/MMP, plots of these 
functions are generally referred to as ambiguity surfaces. 

b. Matched-Mode Processing (MMP) 

MMP [Refs. 4, 49, 50, 51] requires the pressure field in the 


channel to be expanded in terms of normal modes. As will become clearer 
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later, it may be understood intuitively as a transformation of the 
observations from the space of individual hydrophones to the space of modal 
amplitudes. This transformation is known as mode filtering. As before, we 
will assume that an accurate representation of the field (23) requires only a 
finite number of terms M. 

As with MFP, the pressure field is sampled with an N-element 
vertical array of identical elements at depths (2,, z,,...,z,,), located a distancer 
relative to an arbitrary origin. The vector of received pressures (defined in 


(23) and (25)), without additive noise, may be expressed in matrix form as 


p-Zu, (31) 
where 
2ı r,Zı) zi) LEA) 
Z- Z(r.a) Aas | i and 
ZEZ ele êlo Za J 
TR 
uci mv 


Yep JER (sen a — i| (a) (32) 


To 


Z thus contains all information about the receiving array, while u contains 
all information about the location of the emitter relative to the receiver. Both 
Z and u contain information about the channel via the eigenfunctions Z,, 
which are derived from normal mode analysis. An estimate t of u may be 
obtained from (31) using either of two classes of methods. 

The first class of methods regards estimation of à from a purely 


mathematical perspective, namely, as a least-squares problem [Ref. 33], 
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either overdetermined (N2M) or underdetermined (N«M). In this problem, à 


is selected to minimize the quantity 
^ 2 
|Zà - p|;., (33) 


the squared Euclidean norm of the residual (Zà-p). In the overdetermined 
case, assuming that Z has full rank, the solution à is unique. In that case, 
the selection of a suitable method 1s based on considerations of numerical 
stability and computational complexity. If the problem is underdetermined 
(as is the case with the data to be analyzed in this dissertation) or Z is rank- 
deficient, there is an infinitude of vectors à which minimize (33). Two 
subclasses of least-squares methods exist in this case: those which produce a 
solution à with minimum norm (such as the pseudo-inverse method discussed 
below) and those which do not (such as certain versions of the QR 
factorization method). Methods in the latter subclass give solutions with 
significantly greater sensitivity when p 1s contaminated by observation noise. 
We will consider only the pseudo-inverse method in this dissertation. 

The singular value decomposition (SVD) of Z 1s given by [Ref. 
33] 


U?ZV - diag(c,,0,,...,0,) 


0, 


| 
'0l-x 
e 


where U and V are NxN and MxM unitary matrices respectively, 
p=min(M,N), and the g, are the (real and non-negative) singular values (some 


of the o; will be zero if Z is rank-deficient). Note that the matrix partition 
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shown assumes NSM; the partition when N>M is analogous. The 


pseudoinverse Z' of Z 1s defined as 


Zt = VXU”, 
where 
5 = diag dere) (34) 
Oi On 


and R=rank(Z). The minimum-norm solution of (31) can be expressed as 


ü-Z'p. (35) 


Obviously, the expression (30) is very sensitive to the presence of small but 
non-zero singular values (due, e.g., to roundoff error). This phenomenon may 
be satisfactorily dealt with by treating all singular values on the order of 
machine precision (or smaller) as zero. As 1s well known, if N2M and Z has 


full rank, the unique least-squares solution is given by 
a =Z'p =(Z"Z) Z*p. 


This method of modal decomposition is referred to by Yang [Ref. 52] as the 
“Eigenvector Method”. Obviously, modes which are so poorly sampled by the 
receive array that the corresponding columns of Z are nearly linearly 
dependent cannot be resolved using this approach. 

As mentioned earlier, the analysis to follow does not assume the 
existence of an estimate of the noise covariance. It is worth noting that, when 
p is contaminated by observation noise with known covariance (1.e., known, 
except possibly for a constant multiplicative factor), à should be determined 


via the generalized least-squares method [Ref. 33] in order to ensure that 
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undue weight is not given to data from hydrophones with high levels of 
observation noise. Since this method will not be used in our analysis, it will 
not be discussed further. 

The second class of methods (see, for example, [Ref. 51]) for 
estimating à regards the problem from a more physical perspective; i.e., it 
exploits the orthogonality property of the mode functions (columns of Z) (21). 
In (21), if the spacing of hydrophones is sufficiently dense, the integral may 
be replaced by a sum without significant loss of accuracy. Furthermore, 
within the water column, the density p is approximately independent of 
depth. Thus, the columns of Z are approximately orthogonal, at least for the 
lower modes (i.e., for those modes which are well sampled by the receiver 


depths z,), i.e., 
ZEZ = £ diag(v,, Va». Vu), 


where Az is the spacing of the vertical grid on which the mode function is 


evaluated. To obtain an estimate of u, we premultiply (31) by Z” 
Z”p = Z” Zu 


P 
= —diag(V,, V>,...Vy) u. 
A g( 1: V2 m) 
We may thus take the estimate of u to be 


e ED (36) 


This method of estimating u, which we will refer to as the projection method, 
is obviously very similar mathematically to the pseudoinverse method. 


Because of the mode orthogonality assumption, only modal amplitudes 
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corresponding to those columns of Z which form a (nearly) mutually 
orthogonal set may be accurately estimated. 


When observation noise is present, (31) becomes 
p= Zu+n. 


As mentioned earlier, because we do not assume that an estimate for the 
noise covariance is available, the presence of this observation noise does not 
affect how we estimate u (although, of course, the value of the estimate will 


be affected). Using the pseudoinverse method to estimate u, we obtain 


ù = Z'p 

= Z'Zu+Z'n 

=u+Zn . (37) 
From (32) we have 

u = as(t), 
where 
a(roz sa (rz rans MSIE (38) 
Gee TIAE ZA SER) os] i Tn and 
s(t)zexp(j ox). 

If we define 
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(37) becomes 


ù = as(t)+ n’. (39) 


The expression (39) is of the form (1), again for the special case of a single 
signal. Therefore, once à is known, the MFP techniques discussed above may 
be used to find the emitter range and depth. As discussed earlier, no estimate 
of the observation noise will be used in the analysis to follow. Thus, the fact 
that the "noise" n' has a different character from n does not affect our 
analysis. The analysis for the case of the projection method is essentially the 
same as the foregoing and will not be presented separately. Yang [Ref. 52] 
notes that this method of modal decomposition gives a localization estimate 
which is mathematically equivalent to the MFP approach when the Bartlett 
processor and all modes are used. 

Regardless of which method of estimating u is used, care is 
required in selecting which subset (i.e., which components of à and a) of the 


full mode set (obtained from either (35) or (86)) is to be used, because: 


e As the mode number increases, so does the vertical wavenumber 
[Ref. 42]; thus, for higher-order modes, a closer vertical spacing of 
receive hydrophones is necessary for adequate “sampling” 


(analogous to the sampling theorem in time-series analysis); 


e Propagation models generally show greater sensitivity to 
uncertainties in the environmental parameters when predicting 


high-order modes than when predicting low-order modes (see, e.g., 


[Ref. 52]); 
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e Modes which are only weakly present (i.e., which have low modal 
amplitude) at the receiver location can generally not be estimated 


accurately; 


e Only those modes which have most of their energy at depths within 
the physical extent of the receiving array are likely to be estimated 


accurately (see, e.g., [Ref. 52]); 


e When the number of receive hydrophones is less than the number of 
modes supported by the channel (as is the case with the data to be 
analyzed later in this dissertation), the inversion of (31) is an 
underdetermined problem and therefore cannot provide accurate 
values for all modes. 

The first and second considerations usually favor the lower-order modes. This 
generalization is not valid in all situations: for example, incorporation of 
poorly resolved higher-order modes into the estimator can sometimes reduce 
sidelobe levels when the low-resolution Bartlett processor is used [Ref. 52]. 
The third consideration also tends to favor low-order modes, at least when 
the emitter and receiver are widely separated (since higher-order modes are 
attenuated more rapidly). The fourth consideration is relatively easy to 
employ, since the modal structure (1.e., mode shapes) at the receiver location 
is known. The fifth consideration favors modes which are well sampled by the 
receiving array and which are therefore nearly orthogonal, since only nearly 
orthogonal modes are well resolved by mode filtering. Again, these modes are 
generally the low-order ones. In general, an initial localization estimator 
should be constructed based on a mode set selected using the above 
considerations. The peaks of this estimator may be regarded as candidate 


emitter locations. Then, revised estimators (one estimator per candidate 
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source location) may be obtained by using only those modes which are 
expected to be present at the receiver due to sources at these candidate 
locations. Thus, a strategy of iterative improvement may be used to refine the 
estimator. 


The MMP technique may be summarized as follows: 


e Perform quadrature demodulation on the received (real) signals at 


each receive hydrophone to obtain p(t); 


e Use either the pseudo-inverse or projection methods to estimate u(t) 


from p(t) using (31); 


e Generate an estimate of the modal covariance matrix R,,-E[uu"| 


(or the 4th order modal cumulant matrix C,); 


e Using a suitable propagation model and invoking the principle of 
acoustic reciprocity, generate the array manifold vectors a(r,,z,) 


(from (38)) for values of (7,z,) on a suitable grid; 
e Select a subset of the full mode set for use in further processing; 


e Generate the functions in equations (27) (Bartlett), (28) (MVM), 
(29) (MUSIC), or (80) (cumulant MUSIC), as desired. The estimated 
emitter location is the (7,,z,) combination for which the function is 


maximized. 


e As is apparent from our discussion above concerning mode 
selection, the major advantage (at least from the perspective of our 
research) of performing MFP in mode space (i.e., MMP) rather than 


in hydrophone space is that estimation errors due to environmental 


4] 


mismatch may be reduced by using by using only robust modes, 


which are less sensitive to such mismatch. 
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III. EXPERIMENTAL SETUP 


This chapter provides an overview of those aspects of the Barents Sea 
Polar Front Experiment which are applicable to this research. This 
experiment provided the data on which the MFP/MMP algorithms were 
tested. 


A. ENVIRONMENT 


The data used in the analysis to follow was obtained during the 1992 
Barents Sea Polar Front Experiment [Refs. 16, 17, 18, 19]. Most of the details 
of that experiment are not pertinent to our analysis, but may be found in the 
listed references; the pertinent aspects are provided below. 

1. Bathymetry 

Figure 7 shows the bathymetry of the acoustic channel, as well as the 
locations of the source (far left side of the plot) and the receiver (located at 
roughly 34 km range) (to be discussed later). The geometric axis from the 
source to the receiver was almost directly downslope. 

2. Sound Speed Profile 

Figure 7 and Figure 8 illustrate the sound speed field in the channel. 
The curves in Figure 8 were obtained by interpolation of sensor casts made at 
roughly 10 km intervals; the dotted and solid curves correspond to the source 
and receiver locations, respectively. The sound speed field obtained from the 
interpolation was used as input to the BBCM model for all data analysis. The 
higher resolution sound speed field of Figure 7 was obtained by tomographic 
inversion [Ref. 19]; it is provided for illustration only and was not used in our 


analysis. 
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Figure 7: Sound speed field and bathymetry 

The plot uses MATLAB’s “pseudocolor” feature: the gray level at each 
point in the plot maps to a sound speed (in meters per second) per the gray 
level bar at the left side of the plot. The plot clearly shows the front which 
was the primary subject of interest in the Barents Sea Polar Front 
Experiment; this front was nearly perpendicular to the axis between the 
source and receiver. This fact, combined with the fact that sound propagation 
was almost entirely downslope, allows us to make the assumption that no 
horizontal refraction or azimuthal scattering occurred [Ref. 19]. The front 
was observed to move upslope and downslope with a dominant periodicity of 
about two hours and a peak-to-peak amplitude on the order of 4 km. The 
sound speed and density in the bottom were obtained from standard Navy 
databases and were found to be 3200 m/s and 2600 kg/m’, respectively; these 


values were verified by SUS measurements [Ref. 53]. 
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Figure 8: Sound Speed Profiles 


B. SIGNAL AND NOISE 


The transmitted signal consisted of M-sequences with a center 
frequency of 224 Hz. Two sets of M-sequences, separated by about nine hours, 
were transmitted during the experiment. Each set consisted of 30 M- 
sequences, each of about five-second duration, giving a total of about 2.5 
minutes of transmission per set. The unique properties of the M-sequence 
were needed to achieve the goals of the Barents Sea Experiment (i.e., 
accurate travel time determination), but are not relevant to this dissertation 
and will therefore not be discussed here. For the present analysis, the 
received signal was filtered using a Minimum-Variance filter to remove the 
M-sequence properties; this type of filter was selected to ensure that any 
strong interfering signals due to engine noise from the test ship were nulled 


out (see, e.g., [Ref. 54]). We found, not surprisingly, that use of this filter 
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resulted in emitter location estimates superior to those obtained using the 
usual Fast Fourier Transform (FFT) approach. Figure 9 shows the power 
spectral density (relative to the peak) of the unfiltered received signal at a 
particular hydrophone and a particular time. The received signal (at 224 Hz), 
when averaged over all hydrophones, has a SNR of about 10 dB (ratio taken 
over the entire signal bandwidth); the exact value of the SNR is not 


particularly important for our analysis. 
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Figure 9: Power Spectrum 


C. RECEIVER 


The receive array consisted of a vertical string of 16 identical, 
omnidirectional hydrophones with 10 m spacing. The uppermost hydrophone 
was at a depth of 123.8 m. Prior to the collection of data used in this 
dissertation, hydrophones 1, 2, 4, 6-8, and 16 experienced partial failure 


(separation of the two piezoelectric cylinders comprising each hydrophone) 
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which reduced their sensitivity by 6 dB. This sensitivity reduction was 


incorporated into the data analysis. 
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IV. ANALYSIS AND RESULTS 


This chapter provides further details of the analytical technigue, 


including preprocessing methods. It also provides and interprets our main 


results (which take the form of ambiguity surface plots) from application of 


the MFP/MMP estimators to the Barents Sea data for various choices of 


parameters (noise, data length, etc.). 


A. ANALYSIS TECHNIQUE 


1. 


Signal Processing 


The computational procedure may be outlined as follows: 


Construct the density and sound speed fields (functions of range 


and depth) from in situ measurements using bilinear interpolation); 


Using sound speed, density, receiver horizontal location, and 
bathymetry as inputs, and regarding each receive hydrophone as a 
unit emitter, generate (using a suitable propagation model) and 
store the parameters appearing in (31) (a separate set of 


parameters for each hydrophone depth); 


Using a numerical implementation of (31) and invoking the 
principle of acoustic reciprocity, calculate the array manifold on a 


suitable grid of (7,,z,) values; 


Take the raw hydrophone data (an M-sequence) and add the 
desired amount of noise (measured in the same environment during 


a period when no transmissions occurred); 
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e Pass the resulting signal through a 50th order minimum-variance 
filter to eliminate out-of-band noise and to remove the M-sequence 
properties from the signal (i.e. make it into an ordinary narrow- 


band process); 


e Calculate estimated spatial and modal covariance and cumulant 


matrices from the resulting signal; 


e Match these covariance matrices against the replica fields 
calculated above (for MMP, use only the components corresponding 
to the desired modes). 

All computations were done using Matlab 4.2c on a Hewlett-Packard 

735 Workstation. The propagation modeling used the Broadband Coupled 
Mode (BBCM) algorithm [Ref. 19]. In every case, the ambiguity surfaces were 
calculated for a grid with the following specifications, selected to ensure 


peaks would not be missed while keeping the computational load reasonable: 


Parameter Minimum Maximum Increment 
(Emitter) Range 15 km 40 km 40m 
Depth 2m 146 m 2 m 


Table 1: Computation grid 


The emitter range used in the plots is measured with respect to a reference 
different from those shown in Figure 6 and Figure 7: it is measured with 
respect to a point 1675 m downslope from the receiver. The 15-40 km range 
window was selected to allow a fair assessment of our estimators, while 
keeping the amount of computation manageable. Because the water gets 
deeper as one gets closer to the receiver, more modes are required to 
construct the pressure field (the number of modes supported is roughly 


proportional to the water depth for a fixed frequency). The run time required 


o0 


by the BBCM model appears to be proportional to the number of modes 
cubed. 

2. Plotting 

The ambiguity surface plots were generated using Matlab (version 
4.2a) on a Macintosh LCIII and printed on a 600 dot per inch (dpi) HP 
LaserJet 4, using the Matlab “pseudocolor” plot function with 32 gray levels. 
Because a 600 dpi printer is not able to generate a dot screen with adequate 
resolution to display data on a grid as large as that found in Table 1 without 
requiring an excessive amount of space on the page, the ambiguity surface 
data was smoothed before plotting: at each depth grid point, the value plotted 
is the average for three adjacent range grid points. Subjectively speaking, 
little information appeared to be lost as a result of this smoothing. 

The quantitative assessment of the effect of various parameters on 
localization performance presents an interesting problem. On all plots, the 
bright areas correspond to maxima of the functions (27) through (30). No gray 
scale is provided, since, for the MUSIC algorithm, the peak height does not 
necessarily correspond to the power of the received signal. In fact, the values 
corresponding to “white” and “black” differ somewhat from plot to plot. A 
more suitable approach for quantitatively comparing the different plots may 
be based on three major performance measures used in DOA estimation 


literature (see, e.g., [Ref. 23]): 
e Bias in the emitter location estimates; 
e Ability to resolve multiple closely spaced emitters; and 


e Existence of peaks at locations where no emitters exist. 
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In the present case, because the precise location of the emitter is not known, 
the first measure is not applicable in the form stated. Also, only one emitter 
is present, so the second measure is not particularly useful either. The third 
measure is applicable, although not easy to quantify. We have elected instead 


to use the following: 


e Ratio of the height of the correct peak to that of the highest false 
peak (M1); 


e Size (area) of the correct peak relative to the area of the entire grid 


(M2); 


e Ratio of the average height of all false peaks to the height of the 
correct peak (M3). 


Obviously, when M1<1, an incorrect estimate for the emitter location will be 
obtained. The values of these three measures are given in the caption for 


each figure that follows. 


B. OVERVIEW OF RESULTS 


In each section to follow, we consider the effect of the variation of a 
single parameter (e.g., data length, algorithm type, SNR, etc.) while holding 
all other parameters fixed at the following nominal values (for which 


performance is good): 


e Matched-mode method (using projection method of mode filtering 


and modes 1-4); 


e Data length of 1024 points; 
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e Data taken from the beginning of the second set of M-seguence 


transmissions; 


e SNR-10 dB (i.e., no noise beyond that actually observed with the 


signal); 
e MUSIC method of array processing; 
e Second-order statistics; 
e Coupled-mode propagation model; and 


e Hydrophone data corrected for reduced sensitivity of damaged 


phones. 


Each plot to follow has a title at the top containing the most 
significant parameters (MFP/MMP, SNR, data length, and array processing 
algorithm). Additional pertinent information appears either in the caption or 
in the text referring to the plot. In every case, the correct emitter location is 


at approximately 36 km range and 122 m depth. 


C. MATCHED-FIELD VERSUS MATCHED-MODE PROCESSING 


As mentioned earlier, propagation models are generally more sensitive 
to errors in knowledge of the environment when predicting high-order modes 
than low-order modes. We therefore expect MFP to exhibit a higher incidence 
of false peaks than MMP, since MFP uses all available modes (26). Figure 10 
(MFP) and Figure 11 (MMP) illustrate this effect. For purposes of 
illustration, we have added a small white circle at the correct emitter location 
in Figure 10, although we have not done so with subsequent plots, because 


the emitter location is fixed. Although MFP gives a peak at approximately 
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the correct emitter location (M1=0.94), there are 2 other peaks which are 
higher, as well as numerous smaller peaks (M3=0.043). MMP gives the 
highest peak at the correct location (M1=1.7) and has fewer and smaller false 


peaks (M3=0.016). 
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Figure 10: Matched-Field Processing (M1=0.94; M2=0.0077; M3=0.043) 


54 


MM, MUSIC, Modes 1-4, SNR=+10 dB, 1024 pts 





20 


40 


Depth (m) 


60 


80 


100 


120 


140 


1.5 2 2.5 3 3.5 
Range (m) 


x1g" 


Figure 11: Matched-Mode Processing (M1-1.7; M2-0.012; M3-0.016) 


D. BEHAVIOR OF DIFFERENT ARRAY PROCESSING METHODS 


Both Figure 12 (MVM) and Figure 13 (Bartlett) show poor resolution 
compared with MUSIC (Figure 11). Both of these methods produce the 
largest peak at the correct location (M1=1.26 and 1.03, respectively), but 
there are large and numerous false peaks (M3=0.11 and 0.14, respectively). 
In particular, with the Bartlett method, the correct peak is nearly impossible 
to identify visually. The behavior of these three methods when used on this 
data set is thus consistent with that observed in DOA estimation and with 


modeled MFP data (see, e.g., [Refs. 3, 23]). 
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Figure 12: Minimum-Variance Method (M1=1.26; M2=0.033; M3=0.11) 
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Figure 13: Bartlett Method (M1=1.03; M2=0.017; M3=0.14) 
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E. EFFECT OF DATA LENGTH ON PERFORMANCE 


As the length of data increases (assuming temporally stationary data 
and fixed SNR), the accuracy of the estimate of spatial covariance should 
improve. On the other hand, if significant non-stationarity is present, we 
expect that increases in data length past a certain point (depending on SNR) 
may actually degrade the accuracy of this estimate. This is the case with our 
data, as may be seen in Figure 14 through Figure 16, where the ambiguity 
surface for a data length of 512 shows some improvement (M1=2.2, 
M3=0.015) with respect to the surface for a data length of 1024 (Figure 11), 
and the size of the false peaks increases noticeably for data lengths of 2048 
(M1=1.6, M3=0.019) and 4096 (M1=1.4, M38=.020). This behavior is 
presumably due to surface wave effects, which are the only likely source of 
temporal variability over the roughly one-second time scale involved here. 
Obviously, data length can only be reduced so far before the benefits gained 
by avoiding non-stationarity are outweighed by estimation errors arising 
from low information-to-noise ratio. In fact, performance degradation became 
evident at SNR=+10 dB when the data length was reduced to 256 points (not 


shown). 
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Figure 14: Data Length 512 (M1=2.2; M2=0.012; M3=0.015) 
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Figure 15: Data Length 2048 (M1=1.6; M2=0.012; M3=0.019) 
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MM, MUSIC, Modes 1-4, SNR=+10 dB, 4096 pts 
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Figure 16: Data Length 4096 (M1=1.4; M2=0.012; M3=0.020) 


F. EFFECT OF NOISE ON PERFORMANCE 


Figure 17 through Figure 19 illustrate the effect of additive noise. 
These plots exhibit high false peaks compared with Figure 11 (no added 
noise). Below a SNR of —10 dB, the MUSIC method using 2nd order statistics 
no longer gives the largest peak at the actual emitter location (M1=0.72 for 
Figure 18). Interestingly, the MUSIC method with 4th-order statistics 
actually gives poorer results than with 2nd-order statistics at all SNRs 
(Figure 19 shows the 0 dB result). One reason for this somewhat surprising 
result appears to be that the signal turns out to be roughly as close to 
Gaussian as the additive noise, as measured by the difference between its 
4th-order moment and the 4th-order moment of a Gaussian process with the 
same lower-order moments (11). Specifically, let us define a measure G for 


quantifying Gaussianity for the received signal as follows: 
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ez El: 
| M, |; 





(40) 


where |||, indicates the Frobenius norm", C, is defined by (14), and M, is the 


matrix of 4th order moments; i.e., 
M,(i, J) = E| p:P; P;P; |. 


As discussed earlier, if the processes p, are Gaussian, C,=0 and therefore 
G=0. We find that, in this experiment, for noise alone, G=0.13 and for signal 
alone G=0.12 (the use of different norms in (40) does not significantly affect 
these values). Recall that our motivation for using higher-order statistics in 
the first place was based on an expectation that the noise would be Gaussian 
and the signal would not. It appears that this Gaussian property is 
characteristic of the filtered M-sequence signal; for continuous wave (CW) 
data gathered later in the experiment (at which time, unfortunately, 


environmental measurements are not available), G is significantly higher. 


" The Frobenius norm of a matrix A is defined as AI pU > EVE , Where A; are the 
i,j 


components of A. 
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Figure 17: SNR = -10 dB, 2nd order statistics 
(M1=1.1; M2=0.014; M3=0.028) 
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Figure 18: SNR - -11 dB, 2nd order statistics 
(M1=0.72; M2=0.012; M3=0.035) 
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MM, MUSIC, Modes 1-4, SNR=+0 dB, 1024 pts 
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Figure 19: SNR = 0 dB, 4th order statistics 
(M1=0.34; M2=0.008; M3=0.015) 


Å more significant flaw in the cumulant-based MUSIC estimator is 
illustrated in Figure 20 through Figure 23, which show the performance of 
cumulant-based MUSIC MFP versus that of conventional MUSIC MFP when 
applied to a synthetic data set (no mismatch between the actual and predicted 
sound fields) for SNRs of +10 dB and 0 dB. At -10 dB, the cumulant-based 
method (M1=12.3) greatly outperforms the conventional method (M1-2.9). 
However, at 0 dB, the conventional method is superior (M1=2.7 vs. M1=1.6), 
despite the fact that the signal subspace is estimated much more accurately 
for the cumulant method than for the conventional method (as quantified by 
the angle between the signal subspaces at +10 dB and 0 dB). This behavior is 
due to the fact that both the cumulant matrix C, and the array manifold 
vectors (aca’) are defined to be real: since both MFP and MMP rely heavily on 


phase information for accurate localization, use of amplitude information 
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alone causes the estimators to be highly sensitive to noise-induced errors in 


the estimate of C,. 
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Figure 20: Conventional MUSIC, synthetic data 
(M1=2.9, M2=0.0091, M3=0.014) 
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Figure 21: Cumulant MUSIC, synthetic data 
(M1=12.3, M2=0.0021, M3=0.0022) 
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Figure 22: Conventional MUSIC, synthetic data 
(M1=2.7, M2=0.0093, M3=0.015) 


64 


MF, MUSIC, SNR=+0 dB, 1024 pts 





20 


40 


Depth (m) 


60 
80 
100. 
120 


140 ; 
1.5 2— —2-5 3 i 3.5 
Range (m) x 10^ 


Figure 23: Cumulant MUSIC, synthetic data 
(M1=1.6, M2=0.0022, M3=0.015) 


G. BEHAVIOR OF DIFFERENT MODE INVERSION METHODS 


Figure 24 shows the result when the pseudoinverse method is used 
instead of the projection method. As mentioned above, the two methods are 
very similar mathematically and give about the same performance (compare 
with Figure 11). The false peaks are slightly larger and more numerous with 
the pseudoinverse method (M3=0.020 vs. 0.016 for the projection method). 
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Figure 24: Pseudo-inverse mode filter (M1=1.7; M2=0.013; M3=0.020) 


H. EFFECT OF ARRAY SHADING 


Figure 25 shows the effect of not including the required sensitivity 
correction for the failed hydrophones (z.e., no array shading). The higher false 
peaks are apparent (M1=1.1 and M3=0.025, as compared with M1=1.7 and 
M3=0.016 when shading is used). 
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Figure 25: No phone sensitivity correction 
(M1=1.1; M2=0.012; M3=0.025) 


I. TEMPORAL VARIABILITY OF RESULTS 


Figure 26 shows the results from a data segment obtained during the 
first set of M-sequence transmissions. Although there is a peak at the correct 
emitter location, it is not the largest peak (M1=0.46). The change in 
localization performance with respect to that obtained with a data segment 
from the second set of transmissions is probably due to temporal fluctuations 
in the sound speed field, since none of the other physical parameters changed 
significantly between the two data sets. This behavior is not surprising, since 
the period of frontal motion (two hours) is much less than the time between 
the sets of M-sequence transmissions. Although the plots are not shown here, 
we found that the localization performance remained qualitatively the same 


for all data segments within a given set of transmissions. 
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Figure 26: Data from Ist transmission (M1=0.46; M2=0.015; M3=0.026) 


J. EFFECT OF MODE SELECTION ON MMP RESULTS 


As discussed above, it is important to choose a suitable mode subset to 
ensure satisfactory results. Obviously, with 60 modes available, there is a 
very large number of potential combinations, only a few of which will be 
presented here. Figure 27 shows the result when only the first 3 modes are 
used; a small peak is visible at the correct location (M1=0.46). Figure 28 
shows the result when modes 1-5 are used; the plot shows a slight 
improvement in performance compared to the result with modes 1—4 (M1-1.8 
and M3=0.015 vs. 1.7 and 0.016, respectively). Using more than the first five 
modes tends to increase the number and size of the false peaks. Figure 29 


(modes 1-6 used) shows this effect (M1=1.3, M3=0.018). 
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Figure 27: Modes 1-3 (M1=0.46; M2=0.0070; M3=0.025) 
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Figure 28: Modes 1-5 (M1=1.8; M2=0.011; M3=0.015) 


69 


MM, MUSIC, Modes 1-6, SNR=+10 dB, 1024 pts 
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Figure 29: Modes 1-6 (M1=1.3; M2=0.011; M3=0.018) 


K. EFFECT OF MODEL SELECTION ON PERFORMANCE 


Figure 30 shows the result when the mode coupling accounted for by 
the BBCM method is not included in calculation of the replica fields, that is, 
when the range dependence of the replica fields arises only from the range 
dependence of the mode functions and wavenumbers (the adiabatic 
approximation of (22)). There is no indication of a peak at or near the actual 
emitter location, so the performance measures M1, M2, M3 are meaningless 


and are not provided. 
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Figure 30: Adiabatic model 
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V. CONCLUSIONS 


The results presented in the body of this dissertation clearly 
demonstrate that MUSIC-based MMP techniques may be effectively 
employed even in a very challenging acoustic environment such as that which 
existed during the Barents Sea Polar Front Experiment. It is appropriate to 
review some of the unique features of this environment (as compared to 


idealized numerical simulations or simple experiments): 


e Strong range dependence of bathymetry and the sound speed field 
(which included a strong, rapidly moving front); 


e A degraded receive array spanning about half the water column 
and having many fewer elements than the number of modes 


supported by the channel; 


e Relatively coarse sampling of the sound speed field (only about once 


per 10 km interval, notwithstanding the presence of a front); 


Despite these challenges, we achieved considerable success with our 
localization approach. Some of the specific findings and original contributions 


of this research are: 


e Contrary to much conventional wisdom, the subspace-based MUSIC 
method produced good results despite the inaccuracies inherent in 
experimental data. In fact, the MUSIC algorithm’s high resolution 
was vital to accurate localization, because the receive array was 


relatively short and few robust modes were available. 
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e The Broad Band Coupled Mode (BBCM) model generates replica 
fields with sufficient accuracy to allow localization via MMP in this 
strongly range-dependent environment. The adiabatic approxima- 


tion was found to be grossly inadequate for this environment. 


e The cumulant-based MUSIC estimator used in this dissertation 
was too sensitive to noise-induced and model-induced estimation 
errors to be useful with real data. This behavior was due to the fact 
that the cumulant matrix and the replica fields were defined as real 
quantities; thus, the phase information, which is vital to robust and 


accurate localization, was not available. 


e In an environment with strong temporal variability, localization 


performance can vary drastically over relatively short time scales. 


In summary, our approach, which combined the high-resolution MUSIC 
algorithm with MMP, allowed accurate localization estimates, even though 


only a few robust modes could be obtained via mode filtering. 


The approach used in our analysis may be modified in three obvious 
respects, each of which offers potentially significant improvement in 


localization performance and is worthy of further study. 


The first modification relates to the assumptions concerning 
observation noise. The basic MUSIC algorithm used in our research assumes 
that the noise covariance is some multiple of the identity matrix (i.e., 
spatially isotropic). As mentioned earlier, the MUSIC algorithm may be 
extended to the case where the noise is not isotropic via use of a generalized 


eigendecomposition. Use of this extended MUSIC algorithm has the potential 
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to further lower the SNR threshold (the SNR above which the estimator 


performs satisfactorily) observed in our results. 


A second way to improve localization performance involves 
modification of our cumulant-based MUSIC estimator. We noted earlier that 
an extended version of the cumulant matrix has been defined [Ref. 20]. 
Although an estimator based on this matrix would be more computationally 
intensive than the estimator defined here, it is expected to be less sensitive to 


estimation errors. 


Refinement in the process of selecting suitable modes for MMP 
localization is a third means of improving the estimator. Although the simple 
approach described here was effective in generating an appropriate mode set 
for this environment, it may not produce results of the same quality in other 
environments. An approach relying more heavily on the propagation physics 
of the channel could greatly reduce the amount of “trial and error” involved in 


the process. 
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